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1 Introduction and motivation 

Nonlocal models have been recently used in many applications, including continuum 
mechanics [5T, graph theory [26] . nonlocal wave equations (40] . and jump processes 
[4j HOj ; we consider nonlocal diffusion operators which arise in several and diverse 
applications such as image analyses [§1 US [55J [55] , machine learning [35] , nonlocal 
Dirichlet forms [2] , kinetic equations [3J [51] , phase transitions [5j [TB] , nonlocal heat 
conduction 7] , and the peridynamic model for mechanics [121 137j . In this work we 
consider a nonlocal integral operator for anomalous diffusion which has, as special 
cases, the fractional Laplacian and fractional derivative operators that are com- 
monly used to model anomalous diffusion |29| . Physical phenomena exhibiting this 
property cannot be modeled accurately by the usual advection-dispersion equation; 
among others, we mention turbulent flows [111136] and chaotic dynamics of classical 
conservative systems [15] . 

Nonlocal models differ from the classical partial differential equation models in 
the fact that in the latter case interactions between two domains occur only due 
to contact, whereas in the former case interactions can occur at a distance. In 
particular, let f2 cz W 1 denote a bounded, open domain. For it(x) : £1 — > K, define 
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the action of the nonlocal diffusion operator C on the function u(x) as 

£ U (x):=2f ( M (y)- U (x)) 7 (x,y)dy VxeSicR", 

where the volume of SI is non-zero and the kernel 7(x, y) : SI x SI — > K is a non- 
negative symmetric mapping. We are interested in the nonlocal, steady-state dif- 
fusion equation 

— Cu = f on SI 
u = on Six, 

where the equality constraint (extension of a Dirichlet boundary condition for dif- 
ferential problems) acts on an interaction volume Six that is disjoint from SI. 

The numerical solution of fractional differential equations is an open problem 
and it is the object of current research in applications of models of fractional order; 
see [31] for recent work including many citations to the literature. Common tech- 
niques include methods that take advantage of Laplace and Fourier transforms to 
obtain classical solutions [231 13B] and finite difference methods [S7J 12H1 [Ml 13E1 ST] 
used for constructing numerical approximations. Galerkin discretizations and their 
error analysis have been considered in [17[ 119] for the discretization of the steady 
state fractional advection-diffusion equation. 

A goal of this paper is to develop and analyze discretization methods for frac- 
tional Laplacian equations on bounded domains. We do this by exploiting the fact 
that the fractional Laplacian operator (— A) s is a special case of the nonlocal op- 
erator L. In particular, we compare the solution of the nonlocal steady diffusion 
equation with the solution of the fractional Laplacian equation on bounded domains 
and show that solving nonlocal problems is a viable alternative to solving fractional 
differential equations. A main contribution of this work is to show that not only 
(— A) s is a special case of C, but that it is the limit of the nonlocal operator as 
the nonlocal interactions become infinite, provided that sufficient conditions on cer- 
tain kernel functions hold. This fact has important consequences in the treatment 
of problems involving the fractional Laplacian equations on bounded domains. In 
fact, nonlocal problems are a well posed and are a more general formulation of 
fractional differential models; this is useful for both the analysis of fractional dif- 
ferential equations and for developing finite-dimensional discretization schemes. In 
[15j . finite-dimensional approximations of nonlocal problems is discussed, includ- 
ing finite element discretizations; Galerkin formulations are introduced and results 
are proved about their well posedness and about estimates for the approximation 
error and the condition number of finite element matrices. This helps in design- 
ing efficient numerical methods for the solution of nonlocal diffusion problems and, 
as a consequence, of fractional differential problems. Also, being able to quantify 
the discrepancy between nonlocal solutions and solutions of fractional differential 
equations (as we demonstrate in this paper) allows us to determine to what extent 
the former are accurate approximations of the latter. Furthermore, an important 
advantage of approximating the solution of fractional differential problems with 
nonlocal solutions is that in the latter case we do not have to deal with infinite- 
volume constraints as happens for the solution of fractional differential equations 
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on bounded domains where some expedients for constructing finite-dimensional ap- 
proximations have to be introduced. 

We note that the analysis and approximation of fractional Laplacian problems 
for the case f2 = E™, though of interest in many applications, is not a goal of 
this work and is not discussed. Here, our interest is strictly on bounded domain 
problems. 

We also note that a significant advantage of recasting fractional Laplacian prob- 
lems in terms of the nonlocal problems we consider is that, for the first time, such 
problems can be treated on bounded domains for the case of s < 1/2. Indeed, the 
key to this is our introduction of volume constraints as a generalization of boundary 
conditions; the latter are not well defined for s < 1/2 because traces of functions in 
the energy space associated with fractional Laplacian operators are themselves not 
well defined. 

In our analysis, a recently developed nonlocal vector calculus [14] is exploited 
to define a weak formulation of the nonlocal problem and to prove the convergence 
of the nonlocal solution to the solution of the fractional differential problem. In 
S|2j we provide a brief review of those aspects of the nonlocal calculus that are 
useful in the remainder of the paper, introduce the kernel function, and discuss 
its properties. In S|3j we define the fractional Laplacian as a special case of the 
nonlocal operator C and prove the convergence of the nonlocal operator to the 
fractional Laplacian as the nonlocal interactions become infinite. In £|4] we intro- 
duce finite-dimensional discretizations of the nonlocal problem and we study the 
convergence of the approximate nonlocal solutions to the solution of the fractional 
Laplacian equation. We also discuss the choice of nonuniform grids for mitigating 
the high computational costs that occur when the extent of nonlocal interactions 
becomes large. In f|5j we present results of some numerical tests for finite element 
discretizations of one-dimensional problems; by providing qualitative and quanti- 
tative comparisons between approximate nonlocal solutions and exact solutions of 
the fractional Laplacian equation, these results illustrate the theory introduced in 
§i]and|4} 

2 Elements of a nonlocal vector calculus 

In this section, we review relevant aspects of the nonlocal calculus including non- 
local operators, kernel functions, and nonlocal function spaces. Details about the 
nonlocal calculus are found in [14] , 

The action of the nonlocal divergence operator T> : 1" ->Roni/ is defined as 

2>(i/)(x):=f (i/(x,y) + i/(y,x)).a(x,y)tfy for x 6 E n , (la) 
Jr™ 

where i/(x, y), a(x, y) : R" x W 1 —* E n with a antisymmetric, i.e., a(y,x) = 
— a(x, y), are given mappings. This definition is justified in [T3j. The action of 
the nonlocal gradient operator Q : E" x W 1 — > E n on u is defined as 



0(u)(x,y) := (u(y)-u(x))a(x,y) forx,yeK'\ 



(lb) 
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where u(x) : M. n -> 1 is a given mapping. The fact that — Q and T> are adjoint 
operators is shown in [2]; in fact, in [TJ], the operator Q is denoted by —T>*. 
The operator £ : R n -> K is defined by 

-£u(x) := 23(0 • 0u)(x) = 2 f («(y) - u(x)) 7 (x, y) dy 

Jk™ (2) 

for xeK" with j(x, y) := a(x, y) • (0(x, y) • a(x, y)), 

i.e., it is the composition of the nonlocal divergence and gradient operators, where 
0(x, y) = 0(y, x) denotes a second-order, positive definite symmetric tensor. Note 
that the kernel 7(x, y) is symmetric, i.e., 7(x, y) = 7(y,x). The operator — C is 
non-negative because is positive definite and T> and — Q are adjoint operators. If 
is the identity tensor, C can be interpreted as a nonlocal Laplacian operator. 
The interaction domain corresponding to an open subset fl c R™ is defined by 

fix : = {y e R n \fl such that a(x, y) ^ for some x e fl}. 

Thus, fix consists of those points outside of fl that interact with points in fl. The 
situation fix = which we consider here, is allowable as is the case SI = 1" 

which we do not consider due to our interest in fractional Laplacian problems on 
bounded domains. 



2.1 The kernel 

We assume that the domain fi, as well as fix and f2 u Oi, are bounded with piece- 
wise smooth boundary and satisfy the interior cone condition. For the kernel 7, we 
assume that it is symmetric and satisfies 

7 (x,y)^0 VyeBj(x) 
7(x,y)^ 7o >0 VyeB A/2 (x) (3) 
7 (x,y)=0 Vye (f2 u r2 x )\B A (x) 

for all x e £1 u fix, where 70 and A are given positive constants and B\(x) := {y e 
fl u fix '■ |y ~ x l ^ A}; thus, nonlocal interactions are limited to a ball of radius A 
which itself is referred to as the interaction radius. Note that then 

fix = {y e R n \fl : |y — x| < A for x e fl}, (4) 

i.e., the interaction domain fix is a layer of thickness A surrounding fl. We further 
assume that 

7l , ^ 7 (x, y) 72 , ^ for y e B A (x) (5) 

y — x \ n + 2s ' y — x n+2s 

for all x e SI, where s e (0, 1) and 71 and 72 denote positive constants. An example 
is given by 

7(X ' y) = |y-x|"+^ (6) 
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with cr(x, y) bounded from above and below by positive constants. 

In [TS], other choices for the kernel are discussed and analyzed, showing that the 
nonlocal diffusion operator C has more general applicability than to just fractional 
Laplacian problems. 

Note that we allow for all s e (0, 1) and not just s > 1/2. 

2.2 Equivalence of spaces 

We respectively define the nonlocal energy semi-norm, nonlocal energy space, and 
nonlocal volume-constrained energy space by 

1/2 

l\ f G(v)(x,y) • (e(x,y)-0(«)(x,y))4ydx) (7a) 

V(fiufi z ) := {ve L 2 (Q,kjVl i ) : |||«|||<oo} (7b) 
V C {VL u Ox) := {i> e udj) : i> = on . (7c) 

In [15], it is shown that, for kernels satisfying ^ and ([5]), the nonlocal energy 
space V(fi u fix) is equivalent to the fractional-order Sobolev space H s (fl u ^z)Q 
This implies that V c (£l u Ox) is a Hilbert space equipped with the norm ||| ■ |||. In 
particular, we have 

Cill^ll^(Ounx) < HMH < C 2 \\v\\ HHnuni) Vwe^ufij) (8) 

for some positive constants C% and C%- As a consequence, any result obtained 
below involving the energy norm 1 1 1 • 1 1 1 can be immediately reinterpreted as a result 
involving the norm || • \\ H °(nun x )- 

3 Relations between the nonlocal Laplacian and 
the fractional Laplacian 

In this section, we introduce the fractional Laplacian operator as a special case of 
the nonlocal operator C. Then, we show that for the special kernels given by ^ 
with a = constant, the nonlocal operator C approaches the fractional Laplacian as 
the extent of nonlocal interactions increases, i.e., as A —* oo. 

The fractional Laplacian is the pseudo-differential operator with Fourier symbol 
T satisfying [5] 

7X(-A)V)(0 = |£| 2 m 0< SSS 1, (9) 



For s e (0, 1) and for a general domain O e R™, let 

(«(y) - f(x)) 2 



dydx. 



JnJn |y-x|"+ 2s 

Then, the space H S (U) is defined by [T] H S (D.) := yv e L 2 (f2) : IMI^jq) t I"Ijj3(q) 
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where u denotes the Fourier transform of u. Using Fourier transforms, it can be 
shown [2] that an equivalent characterization of the fractional Laplacian is given 

by 

{ - Ayu= ^ s L T-x\^ dy ' 0<s<1 ' (10) 

where the normalizing constant c n . s is defined as 

Y ( "+ 2 ) 

c "' s = s22s r(i)r(i- s )- 

Here, r(-) is defined via the improper integral 

r(a) = t^e^dt 
Jo 

for all complex numbers a having positive real part. Thus, if the kernel 7(x, y) in 
([2]) is defined as 

7(x ' y)= 2|y-xr^ Vx,yeM«, (12) 

then 

-£=(-A) s , 0<s<l, (13) 

establishing that the fractional Laplacian is a special case of the operator C for the 
choice of 7(x, y) proportional to l/|y — x| n+2s . 
We consider the volume constrained problem 

{— Cu = / in O 
u = in R n \n, ^ ' 

where L is defined as in ^ with 7(x,y) given by (12 1 and / e L 2 (fl). A weak 
formulation of ( 14 ) is 

°f \ Rn £ y| r|" } (v{y) «(x)) dy rfx = £ / 1, rfx V«eff5(R»), (15) 

where H^(R n ) := e H s (R n ) : w = in R"\fi}. The existence and uniqueness 
of such u e H^(R n ) for all s £ (0, 1) is proven in [35]. 
Next, for all x e W 1 , we define 

*(x.v)= \ 2|x-y|« + ^ y£i?A(x) (16) 
yeM»\B A (x). 

Letting £ denote the nonlocal operator associated with 7, we consider the problem 

(17) 

h = in fix, 
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where, again, the interaction domain Six is defined as Six = {ye R"\S1 : |x — y| 
A Vx e SI}. A weak formulation of (17) is 



J I 



2 JouOx J(Slufii)nB A (x) l x — Yl 



2(y) — S(x) ... , , 



\ fvdx V v e V c (n u Sl x ) 



(18) 



where V^(0 u Six), isomorphic to i?Q(K n ) [33], is the volume constrained energy 
space associated with the kernel 7. 

We now show that the solution u of ( 15 ) is the limit, as A —* go, of the solution 
u of (pi). 



Theorem 3.1 Let u e H^(W l ) and u e V C (S1 u Six) denote the solutions of (15) 
and (18), respectively. Then, with I := min{i? : SI c _B fi (x) Vx 6 SI}, 



K„ 



K„ 



(19) 



||t»-«||ia(n) < C 2 s(A _ /)2s HU 2 (f2) 



where C\ is the equivalence constant in (|8| and K n , independent of s and A, is 
determined by the spatial dimension n. 



Proof. Combining ([15]) and ( 18 1 we have 
u(y) - u(x) 



r ^ 

« Jr" |x - 



|n+2s 



(w(y) - u(x)) dydx- 



f f ,rjl ) My)-<x))dydx = V^e^SluSlx) 

J^uOx J(nun x )nB A (x) l x — y| 

If we define Sl c = ]R"\(S1 u Six), this equality is equivalent to 

I J 

\ I 

JQuSlx J(flu!)i)nBj(x 

J I 

I J 

-[ f 



«(y) 




u(x) 


|x- 


y 


n+2s 


«(y) 




u(x) 


|x- 


y 


n+2s 


«(y) 




u(x) 


|x- 


y 


n+2s 


w(y) 




u(x) 


|x- 


y 


?i+2s 


a (y) 




u(x) 



My) 


— v(x)) dy dx+ 


(«(y) 


— v(x)) dy dx+ 


(«(y) 


— u(x)) dy dx+ 


(«(y) 


— u(x)) dy dx+ 


My) - 


- w(x)) dy dx = 
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LI w-yi-Jg-™^-™*** 



1 X J(OuO x )nB A (x) 



I I 



(nuf2z)\s A (x) |x-y 

u(y) - u(x) 



u(y) — u(x) ... , ,, , , 

_ v m-U ; Ky)-«W)dyrfx+ 



LuflxLe |x-y|"+ 2s 

f f u(y) - 

Jo, Jl 



(u(y) - w(x)) dydx+ 



Thus, 



I I 



OuO x J(OuO x )nB A (x) 



|x- y\ n + 2s 
(u(y) - u(y)) - Q(x) - S(x)) 

| x _ y |n+2s 

u(y) - u(x) 



I J 



f f ^ 

|Jnun x Jn c l x 

II I 



|x — 


y 


n+2s 


u(y) 




u(x) 


|x- 


y 


n+2s 


«(y) 




u(x) 


|x- 


y 


n+2s 



(t)(y) - u(x)) dydx = 0. 

(u(y) - d(x)) dydx 

(u(y) - u(x)) dydx -+ 
(v(y) - v(x)) dydx 



(20) 



(v(y) -v(x)) dydx 

We treat Xi, X 2 , and X 3 separately. We first split X\ in the two integrals X a and X,: 



J I 

JouO x J(( 

II 

Jo J(f 

/ I 

Jo x J(f 



u(y) 


- w(x) 


|x- 


y |n+2s 


«(y) 


-u(x) 


|x- 


y \n+2s 


w (y) 


-u{x) 


|x- 


y|n+2s 



(v(y) - v(x)) dydx 
(v(y) - v(x)) dydx 
(v(y) -v(x)) dydx 



= X a + lb- 



Then, 



u(x)v(x) 
Jo J( 

u(x)v(x) 
Jo Jb a+/ 



1 



(OuOx)\Ba(x) |x-y|"+ 2s 



(x)\b a( x) |x-y|"+ 2 * 



dydx 



dy dx 



The value of the inner integral, say i, depends on the dimension n of the problem: 



n = 1, i 



sA 2s 



n = 2, i 



2tt 

sA 2 * 



n = 3, i = 



47T 

sA 2 * 



(21) 
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With k n = 2, 2tt, 4tt for n = 1, 2, 3 respectively, we thus have 



sA 2s 



w(x)i;(x)dx s£ -r?r||w||ia(n)||«||i,a(n) 



(22) 



where we used the Cauchy-Schwarz inequality and ([8|. Next, 

u(y)«(y) 



/iji j(f!ufii)\Bj(x) i x — yi 

w(y)«(y) 



J I 

JB A+J (x)\B x __f(x) Jf 

"(y)^(y) 

Jn Jb x+i 



n |x- y 



n+2s 



dy dx 
dy dx 



1 



+J (£)\s A _ J (a) |x-y|"+ 2s 



dx dy 



where x can be any point in f2; this implies that we can choose x = y for any y e Q. 
Thus, we have 



u(y)v(y) 
Jn Jb x+i 

u(y)v(y)dy 

Jn 



s(A - 1) 2 " 



1 

(y)\B A _ f (y) l x -yl" +2s 



dx dy 



(23) 



Cis(A - J) 



2.s 



IMlL 2 (n)IIMII> 



where the value of the inner integral is obtained in a way similar to (21 ). Next, we 
treat I2 and Z3. 



dy dx 



J u(x)u(x) 

r it 

J ii(x)v(x)dx < Cig p s HU'(n)llfo 



(24) 



Again, the value of the inner integral is obtained in a way similar to (21 1. The 
same bound can be found for I3. Now, letting v = u — u and combining (20), (22 1, 

4fc„ 



((231, and pi}, we have 



u — it < 



Cis(A - /) 
Letting K n = 4fc„ we conclude that 



it — u sS 



K, 



c iS (\-i)^ Ml2 ^- 



The bound for the H s norm is obtained using the equivalence relation ([8]); the 
same bound holds also for the L 2 norm because ||ie|U 2 (n) sg ||ui|liP(OuSiz) for all 
we H s (Quft x ). □ 
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According to this theorem, lower values of s lead to slower convergence, i.e., the 
effects of nonlocality become more pronounced in the sense that the value of the 
interaction radius A has to be bigger in order to achieve the same error for smaller 
values of s compared to larger ones. 

It is tempting to think that the L 2 error estimate in ( 19 ) is pessimistic and 
that a refined analysis would show that that error has a better rate of convergence 
compared to the H s error. However, the numerical illustrations of { 19 show that 
the L 2 error estimate in ( 19 1 is indeed sharp. 



4 Finite-dimensional approximations 

In this section, we consider the convergence of solutions of finite-dimensional dis- 
cretizations (including finite element approximations) of the nonlocal problem to 
solutions of the fractional Laplacian equation. We also discuss the use of nonuni- 
form grids for mitigating the computational costs caused by the increase in the size 
of the interaction domain as A — > oo. 

We consider conforming finite-dimensional subspaces 

V N (Cl u fi z ) c V(tt u fl x ) (25) 

parametrized by an integer N — > oo and then define the constrained finite-dimensional 
subspace 

V C N (Q u Qi) := V N (n u n x ) n V c (fi u fi z ) 

= {v e V N (fl uflj) : v = on fl x } . ' _<> ' 

The common choice for N is the dimension of the subspaces V N (fl u f2j). We 
assume that, for any function v e V(f2 u £l x ), the sequence of best approximations 
with respect to the energy norm ||| • ||| converges, i.e., 

lim inf |||u-«jv||| =0 VveWJ]uJ]i). (27) 

We also define the Galerkin approximation ujy e V^(Cl u Vt x ) of the solution u of 
(18) as the solution of 

z Jnufiz J(nuO x )n-B A (x) |y — x l Jn (28) 

Vu e v^cn u n x ). 

In [15] it is shown that um — » S for oo for all / e L 2 (f2) and s e (0, 1). 



4.1 Finite element approximations 

We consider finite element approximations for the case that both VL u Ox and fi 
are polyhedral domains. We partition O u £l x into finite elements and denote by h 
the diameter of the largest element in the partition. We assume that the interface 
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f2 n fix consists of finite clement faces and that the partition is shape-regular and 
quasi-uniform [5] as the grid size h — > 0. We choose the subspace V^(Cl u f2x) cz 
V^(0 u fix) to consist of piecewise polynomials of degree no more than m defined 
with respect to the partition. We have that N — » oo as h — * 0, where N denotes 
the dimension of the approximating space V c , (fi u fix)- 

In [TS] it is shown that, assuming that u e H m+t (p, u fix) and s ^ t ^ 1, for 
the kernel given in ([6]) and for sufficiently small h, there exists a constant C3 such 
that 

\\u - SAr|| ffs(nu0l) < C 3 /i m+t " s ||S|| /f m+ t(nun:r ). (29) 

Our goal is to find an estimate of the error — 2/v,A||ip(nuSiz)j i-C, of the difference 
between the solution of the fractional Laplacian equation and the finite element 



solution of ( 28 ) for a particular value of the interaction radius A. The estimate is 



easily found combining (191 and (29). By the triangle inequality, we have 



||u - UN,\\\H*(nun x ) < \\u - u\\h'(QuQ x ) + \\u - Sjv,A|ir*(nun x ) 

< a(A ?n2. l"IU»(o) + c s h m+t -*\\u\\ Hm+HnuQx) , (30) 

where C4 = K n jC\. Note that the convergence to the solution of the fractional 
Laplacian equation depends on both A and N; for very fine grids, i.e., very large N, 
the finite element approximation error is negligible whereas, for very large interac- 
tion domains, i.e., very large A, the finite element approximation error dominates. 

Note also that this result holds for a uniform grid on f2 u fix, but, in practice, 
for very fine uniform grids and very large interaction radii, computing Sjv,a is not 



affordable. However, the contributions of the integrals in ( 28 ) over regions of Qx 
far from the domain f2 are not as significant as those of regions close to it; for this 
reason, it might be convenient to utilize a coarser grid as we move further away 
from f2. 



4.2 Nonuniform grid in Vtj 

Let the partition of f2 u fix be uniform in O and nonuniform in f2x and let problem 
( [28| be further discretized by approximating the integrals via numerical integra- 
tion^] we denote by the corresponding solution and find a bound for the error 

\\u - Siv||jj»(nunx)- 

Let A(-, '):Vxl/->M and F(-) : V — > R respectively denote the bilinear form 



and functional associated with problem ( 18 ) 



A( W ,v)= C -^\ f W [ y) jW (v(y)-v(x))dydx (31a) 

z JQuQx J(nufii)nBj(x) W ~ x l 

F(v) = \ fvdx for w.v e V. (31b) 
Jo 

2 Problem \28\ has also to be solved using numerical integrations for the case of uniform grids, in 
contrast to case of the Laplacian operator for which the stiffness matrix can be evaluated exactly. 
However, for the nonlocal case and for uniform grids, the integration error can easily be made 
negligible whereas nonuniform coarser grids in f2i may affect the accuracy of the approximation. 
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Let An(-,-) : V N x V N — > M denote an approximation of A obtained by numerical 
integration of (31a) for w, v e V N and let un denote the solution oQ 

A N (u N ,v N ) = F(v N ) \/v N eV N . (32) 

We assume that An is a coercive (with coercivity constant Cg) and continuous 
bilinear form on V N , then, problem ( |32| ) is well-posed; furthermore, we can apply 
the Strang lemma for generalized Galerkin problems to obtain an a priori error 
estimate for the energy norm [TS]. We have 



| u — un 1 1 1 < inf \ I 1 + — I 1 1 1 u — vn I 



v n bV n IV C 5 



J_ , \A(v N ,w N ) - A N (v N ,w N )\ 

+ sup ... ... 

c 5 u, Jv ey JV \{0} HPivlll 



(33) 



Next, we find a similar result for the H s norm. Let vn be any function in 
V N (fl u fix), we first find a bound for the difference \\un — wjv||j?«(nufJi)- We 
have 

C 5 |||wAr - «jv||| 2 *S Ajv(Sjv - v N ,u N - v N ) 

= 4y(ttjv, U N - V N ) - A N (v N ,U N - Vn) 

= F(2jv - v N ) - A N (v N ,u N - v N ) 
= A(u, u N - v N ) - A n (v n ,un - v N ) 

= A(u - v N , un ~ vn) + M v n,u n - vn) - A N (v N , u N - v N ) 

ujvHI \\\u N - v N \\\ + \A(v N ,u N - v N ) - A N (v N ,u N - v N ) 

Letting AA(vn, vjn) = \A(vn, un~ v n) — An(vn, un— v n)\, the previous inequality 
implies that 

in- in ^ 1 An- in , AA(v N ,u N - v N ) 

\\\Un — vn\\\ < 77- \\\u — vn \\\ + 



C 5 V \\\un - vn 

Using the equivalence of norms (|8| we have 



Ci||Sjv - v N \\H°(nun x ) -? |||«jv - "jv| 

C 2 I 1 AA(v N ,w N ) 
-pr\\ u - v N\\H-(nuQ x ) + 7=r SU P — m m — • 

^5 ^5 w N eV N \{0} Ill-will 



Thus, 



n~ H ^ C 2 ,|~ I, 1 AA{v N ,w N ) 

\\1J-N -VN\\H"(n^Q. x ) < -FT \\ u ~ v N\\H-{aun x ) + FTFT SU P m Ml ' 

°1°5 WU 5 u, jve y«\{0} 1 1 |tf JV 1 1 1 

3 Of course, in practice, one also has to use numerical integration to approximate the right-hand 
side functional F(v) . However, this term does not pose any problems so that we assume it can be 
integrated as accurately as one needs for it not to affect the accuracy of approximate solutions. 
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Now, by the triangle inequality, for all e V N (Sl u fii), we have 
||2 _ "jv!l.ff 3 (nufiz) ^ II" _ v n\\h s (QuQx) + II - u Jv||j?«(nuni) 



sup 

5 w n gV n \{0} 



AA(v N ,w N ) 



\w N \ 



Next, we consider a finite element discretization; we denote by {Ei}^ , {Ei)f =1 a 
n and c the elements of the partition, and by hi the diameter of 

each element; we also denote the (uniform) grid size in fi by h = hi, for i — 1, . . . N. 
Note that 



N+K 



N+K 



< 



2 \\(u-v N )\\ Hs(Ei) . (34) 



i=l 



We conclude that 

/ C 1 \ N+K 

\\u-U N \\ H s (n ^ ax) < ^1+ (T^rJ 2 ll( M -«AfJllif s (B s ) 

AA(t)Af, Wat) 



+ 



sup 



C1C5 u, Ne y«\{o} III^JvlH 
Because the previous inequality holds for all v n e V N , we can write 

I" - u N \\ H s^ nx) < 

N + K 1 A A I \ 

2j IK u_ v at)||jj S (e 4 ) + 

1 



< inf 

v N eV N 



< 11 + 



< 1 + 




sup 



ff s (Ei) + 



inf 



C1C5 u. jv ey«\{0} III^Jvlll 



sup 



C1C5 «»ey tts6 v ff \{0} lll^lll 
1 - .. Aifvjv.iojv) 



inf 



u\\H>m + ^r r — 



sup 



(35) 



where we used (29) and the fact that u\ 



for all 



N + 1, . , . N + K. 



Note that the first term of the right hand side of ( 35 ) depends on the grid size 
inside of fi; thus, it is not affected by a coarser nonuniform grid in fix- This is 
not the case for the second term which depends on the integration method and on 
the grid size inftu fix- In the limit cases of very accurate integration formulas or 
of a fine grid in fix, the second term in (35) is negligible, whereas as the degree 
of accuracy of the numerical integration becomes smaller or as the grid becomes 
coarser, the leading term in the estimate is the "error" AA. Another important 
factor to be considered is the fact that the integrand vanishes in regions of fix very 
far from fi; thus, for such regions a coarser grid might not affect the accuracy of 
the solution. Hence, in order to find a bound for AA (and, as a consequence, a 
coarsening rule) one has to take into account the interaction of the accuracy of 
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the quadrature rule, the grid size, and the structure of the integrand in fix- As 
of now, we do not have such a bound nor an optimal coarsening algorithm, but 
we empirically choose larger hi's far from f2. In §5.1[ we present an example of a 
coarsening algorithm which preserves accuracy, provided that a proper empirical 
tuning of its parameters is performed. 



5 Numerical tests 

In this section we present the results of numerical tests for finite element discretiza- 
tions of one-dimensional problems. Although preliminary in nature, these results 
allow us to illustrate the theoretical results presented in §Sj3] and [4] and to demon- 
strate the viability of using nonlocal diffusion problems as a vehicle for determining 
approximate solutions of fractional derivative problems posed on bounded domains. 
We first consider the convergence of approximate nonlocal solutions to the corre- 
sponding nonlocal solutions and we discuss the choice of a coarser grid in fix in 
terms of accuracy and memory requirements. Then, we provide qualitative and 
quantitative comparisons between approximate nonlocal solutions and exact solu- 
tions of the fractional Laplacian equation. 



5.1 Finite element approximation of a one-dimensional prob- 
lem 

We define a one-dimensional model problem and introduce its finite element ap- 
proximation. Let f2 = (—1, 1) and Six = ( — 1 — A, —1) u (1, 1 + A) so that f2 u fix = 
(— 1 — A, 1 + A). In this case (28) becomes 



l + A rX + X 



2 J-i-aL 



un(v) - u N (x) 



-x \y-x\ 1+2 



(y(y) - v(x)) dydx = J 



fvdx, VveVf. (36) 



For integers K, N > 0, we introduce a partition of f2 u f2x = [—1 — A, 1 + A] such 
that 

- 1 - A = x- K < ■ ■ ■ < x^i < = x < xi < ■ ■ ■ < x N ^i 
< xn = 1 < xn + i <■■■.< xn + k = 1 + A. 

Then, we define h as max i= _x....,A'+A r -i l^i+i — a?* | - We choose V c r (f2 u fix) to be 
the finite element space of continuous piecewise-linear polynomials with respect to 
the partition (37). As a basis, we choose the standard "hat" functions {<pj{x)}^ 



K+N 



-K ' 



Let un(x) = ~Y^=^KUj4>j(x) for some coefficients Uj] because the hat functions 
satisfy 4>i( x j) = we have that Uj = u^j(xj). Then, (36) is equivalent to 



t % u 4Z£-l X y) -J^ ] {Mv) Ux)) dydx = f i M 



(38) 



l,...,N-l, 



Also, because u = e f2x, U$ = Un = 0. 
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The assembly of the matrix and right-hand side of the linear system correspond- 
ing to ( 38 ) requires the use of quadrature rules for the approximate evaluation of 
the integrals involved. We first break up all integrals into sums of integrals over the 
subintervals (xi, Xj+i) of the partition (37). For the single integrals over fl = (—1,1) 
in ( 38 ) , we use a four-point Gauss quadrature rule with respect to each subinterval. 
The double integrals must be treated with greater care, due to the singularity of the 
integrand. For the inner integrals, we split the integral over (x — A, x + A) into the 
sum of integrals over (x — A, x) and (x, x + A), then split each of those integrals into 
a sum over the subintervals of the grid, and then use a four-point Gauss quadrature 
rule over each subinterval; in this way, we avoid the singularity at x. An accurate 
and efficient choice of quadrature rule for the outer integral is a more delicate issue. 
For this purpose, we use the built-in Matlab function quadgk.m which implements a 
high-order global adaptive Gauss-Kronrod quadrature formula which is well-suited 
for integrands having mild end-point singularities |34j . Another advantage of using 
this Matlab function is that all subintervals are processed at the same time so that 
the number of function evaluations is reduced. 

Note that the double integrals on the left-hand side correspond to a nonlocal 
finite element stiffness matrix; an important issue in the numerical solution of 
nonlocal problems is related to its structure. In fact, differently from the local case 
where the matrix is sparse and banded, in the nonlocal case the matrix is dense 
and the inner integral over (x — A, x + A) spans a large part of Q u f2x ; especially 
when A » /, which is often the case. Thus, computations are much more onerous 
compared to that for differential equations. 



5.2 The analytic solution of the fractional Laplacian 



We consider the analytic solution of problem (14) in the ball Br(0) of radius R in 
E". For every function / e L 2 (Br(0)) there exists a unique function u e H s (R n ) : 



given explicitly using the Green's function [T3] by 

«(x) = f G(x,y)/(y)dy, 

JB R (0) 



i(0) 

where 



rn>(x,y) s-i 

G(x,y) = G„, s |x-y| 2 — I (^T^*"' (40) 



with 



(39) 



G -s->° F(g) and rnfxv)- (jR2 - |x|2)( * 2 - |y|2) (41) 

In the particular case of / = 1 in i?i(x) 

r ( -) 

tt(x) = T 2s . +9 ^4 (1 - |x| 2 ) s . (42) 

V ) r (n±2£)r(l + S 1 1U V ' 
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In the one-dimensional numerical tests we consider the following data sets: 



la 



Ila 




lb 



lib 




f(x) = x 
s = 0.40. 



(43a) 



(43b) 



Note that we choose values of s both less than and greater than 1/2. The corre- 
sponding analytic solutions are shown in Figure [l] 



/ = 1 ^ — 
















/ 14, S 


= 0.4 \ 

= 0.75 




Figure 1: Analytic solution u for cases Ia,b in (43a) (left) and cases lla,b in (43b) 
(right). 



5.3 Convergence of finite element approximations to the non- 
local solution 

We examine the convergence, with respect to the grid size h, of finite element 
approximations. These results are not related to the solution of the fractional 
Laplacian equation and they only serve to illustrate the theoretical results about 



the convergence of the approximate nonlocal solutions un to the solution u of ( 18 ). 
We consider the data sets la and lb; because an exact solution is not available, we 
define a surrogate ur for u to be the finite element solution on a very fine uniform 
grid. We then report on the "error" Wur — unWl 2 ^): i-e., the difference between the 
surrogate and the approximation un obtained using coarser uniform grids and the 
same interaction radius A. In particular, the surrogate is determined usinj^A = 0.1 
and h = 2 -11 . In Table ([!]), we list the error and the corresponding rate. In both 
cases we observe a (optimal) rate of of convergence of approximately 0.5. This is 
expected; the analytic solution belongs to H 5 ~ £ (R); in fact, it is continuous but 

4 This small interaction radius has been chosen reduce the computational costs compared to 
that which would be incurred for larger values of A; in fact, in general, with uniform grids on 
O u Qxj one can afford to use only small interaction radii. 
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its derivative has an infinite discontinuity at x = —1,1. 



estimate ( 29 ) does not hold because u £ H m+t (R) with m 
cannot expect a rate better than 0.5. 



As a consequence, the 
= 1 and t e [s, 1] and we 





la 


lb 


h 


error 


rate 


error 


rate 




6.92e-02 




6.23e-02 




2- 4 


4.74e-02 


0.55 


4.53e-02 


0.46 


2~ b 


3.17c-02 


0.58 


3.08e-02 


0.55 


2 -6 


2.11c-02 


0.58 


2.08e-02 


0.55 



Table 1: Dependence on the grid size h of the "error" — un\\l 2 (Q) an d the rate 
of convergence of approximate nonlocal solutions. Results are provided for the test 
cases la and lb. 



Nonuniform grid in Slj For large values of A, using a uniform grid in fix is 
not computationally affordable so that using a coarser grid is mandatory. However, 
this might affect the accuracy of the solution; thus, we need a coarsening algorithm 
that preserves accuracy. Of course, this is not a trivial task (see ( 4.2 1. We design 



an algorithm which is not optimal but still allows us to maintain a certain level of 
accuracy with a proper tuning of parameters. 

Referring to the partition (37), we let x = (xq + xn)/2 and for i 
and i = N + 1, . . . N + K , we define xi as 



-K,. 



1 



Xi = < 



Xi+i - h 



Xi_i + h 



x ~ x i+l 
X — Xq 

Xi— 1 X 

xn — x 



K 



i = N + 1,...N + K, 



(44) 



where K is chosen so that x-k + \ > — 1— A and x^+k-i < 1+A; then, X—k — — 1— A 
and xn+k = 1 + A; h is the grid size inside of £1 and p is the parameter that 
determines the coarseness: the larger p, the coarser the grid. For some values of 
N and A, we test the convergence of the solutions obtained with a nonuniform grid 
to that obtained with a uniform one for several values of p. In Table [2j we report 
the dimension of the partition for N = 2 7 and A = 2 2 , 2 3 ; p = corresponds to the 
uniform grid and it is reported as a reference value. Unfortunately, the difference 
between the solutions on nonuniform and uniform grids, though decreasing for 
decreasing values of p, does not have a specific convergence rate. For the numerical 
tests presented in the following sections, p is chosen empirically in such a way that 
the coarsening does not affect the accuracy, i.e., \AA(w n ,v 
negligible for all wn,vn e V N . Figure [2] shows, for A 
computational grids for p = 0.5, 1, 1.5. 



l/llkivlll 

2 3 and N 



in (35) is 
= 2Vthe 
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p 


X = 2 2 


A = 2 4 


1.000 


337 


413 


0.500 


447 


643 


0.250 


531 


847 


0.125 


583 


985 


0.000 


641 


1153 



Table 2: Dependence on the parameter p of the dimension of the grid for N = 2 7 . 
In this case, for p < 0.5 the coarsening error is negligible. 



5.4 Comparison of u and u 

Through several graphical examples we show how accurately the nonlocal solution 
can approximate the solution of the fractional Laplacian equation posed on bounded 
domains. In Figure [3j left column, we report the finite element solutions un.x and 
the analytic solution u for the data sets la and Ha using several values of the inter- 
action radius A and of the grid dimension N. To appreciate the differences of the 
solutions we report, in the right column, zoom-ins near the peaks of the analytic 
solutions. For a fixed value of N and increasing values of A, the approximate non- 
local solutions converge to u as predicted by theory; we also observe a convergence 
to u for increasing values of N. For all s > 0.5, similar results are obtained, i.e., 
the larger A and N, the closer is un.x to u. 

For s < 0.5, the situation is quite different and deserves a more detailed analysis. 
In Figure[I| we report u N ,\ and u for A = 2 9 , 2 10 , 2 11 and N = 2 7 , 2 8 , 2 9 ; again, the 
plots are zoom-ins around the peaks of u for the data sets lb (top) and lib (bottom). 
Here, changes in A and N have different effects: increasing values of A lead to an 
amplitude reduction whereas increasing values of N lead to an amplitude increase. 
We conjecture that for a fixed very small h, it is possible to observe a convergence 
to u from above for increasing values of A. Because of computational restrictions 
we cannot utilize a grid which is fine enough to observe such behavior. For this 
reason we find it significant to study the convergence of u^.x to u for increasing 
values of A and N simultaneously; in Figure [5] we then observe a convergence to u 
from above. 



5.5 Rate of convergence to the solution of the fractional 
Laplacian equation 

We study the convergence of the approximate nonlocal solution un,x to the analytic 
solution u of the fractional Laplacian equation. Recall that the error bound depends 



on both A and N (see (30)); through several examples, we show how the choice of A 
and N affects the accuracy of the numerical solution and, as a consequence, how it 
affects the rate of convergence. First, we consider simultaneously increasing A and 
N; in Table [3j we report the approximation error, i.e., \\u— UN,x\\L 2 (n), an d its rate. 
In this experiment p = 0.5, except for h = 2~ 8 where p = 0.6. We observe a rate of 
0.5, the same as the one obtained for the finite clement approximations of nonlocal 
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-9 



-1 1 

X 



m 9 • • 



• • • : • 



-9 



Figure 2: Three grid configurations for A 
respectively, from top to bottom. 



2 3 , N = 2 4 , and p = 0.5,1,1.5, 



solutions. We conclude that the convergence with respect to u is dominated by the 



finite element approximation error in ( 30 ) 





la 


lb 


h 


A 


error 


rate 


error 


rate 


2- 4 


2 3 


2.82e-02 




4.47e-02 




2- b 


2 4 


2.00e-03 


0.50 


2.98e-02 


0.59 


2- y 


2 b 


1.41e-02 


0.50 


2.03e-02 


0.55 


2- v 


2« 


9.97e-03 


0.50 


1.40e-02 


0.53 


2 -8 


2 V 


7.05e-03 


0.50 


9.76e-03 


0.52 



Table 3: Dependence on the grid size h and on the interaction radius A of the error 
||u — 2/v,A||z, 2 (n) an d the rate of convergence of approximate nonlocal solutions. 
Results are provided for the test cases la and lb. 



In order to observe the convergence with respect to A, we consider a fixed grid 
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u N (2 4 ) : N = 


2 7 


— u N (2 r '),N = 


2 7 


u N (2 e ), N = 


2 7 


+ u N (2 7 ),N = 


2 7 


— S,v(2 7 ), N = 


2 8 


^S,v(2 7 ), N = 


2 9 


-*-un(2 7 ), N = 


2 io 








5n(2 5 ), N = 2 7 

5n(2 6 ), AT = 2 7 

5n(2 7 ), N= 2 7 
-*-5 w (2 8 ), N = 2 7 
-—U N (2 S ), N = 2 a 
— -5 W (2 8 ), N = 2 e 
-»-ii«(2 8 ). AT = 2 10 




-0.06 


-0.04 -0.02 


0.02 0.04 0.06 






V 








// 










Iff J* 










Iff* 

















Figure 3: Numerical solutions un,\ for different values of N and A show the con- 
vergence to the solution u of the fractional Laplacian as N — > oo and A —* go. On 
the right, zoom-ins of the solutions peaks. Results are provided for the test cases 
la (top) and Ha (bottom). 



and we compare un.\ with it; as we conjectured in the previous section, though we 
do observe a convergence, a (not affordable) very fine grid is required to observe 
the —2s rate. For this reason, we do not find it significant to report those results. 
We examine, instead, the convergence, for a fixed N, of the nonlocal approximate 
solution Un,x to a surrogate ur defined by solving for a finite element approximation 
using a very fine grid and a very large A. In Tables [2] and [5] we report, for the test 
cases I and II, the "errors" |||e||| = |||ur — Ujv,a||| and ||e|| L 2 (n ) = \\u R -u N! x\\L 2 (n), 
i.e., the energy and L 2 norms of the difference between the surrogate and the finite 
element approximation using the same grid and smaller radius. Specifically, the 
surrogate ur for the analytic solution is determined using h = 2~ 9 , A = 2 11 , and 
p = 1.25. We observe the rate —2s predicted by the estimate (19), for both the 
energy norm and the L 2 norm. Note that because the observed rates of convergence 
with respect to the two norms are the same, the sharpness of the L 2 error estimate 
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+ u N (2») 


N = 


= 2 7 


+ un(2") 


N = 


= 2 8 


+ u N (2°) 


N = 


= 2 9 


+ u N (2 10 


, JV 


= 2 7 


+ u N (2 w 


, w 


= 2 8 


+ u N (2 w 




= 2° 


-^E w (2 11 


, w 


= 2 7 


-^S W (2 U 




= 2 8 




, JV 


= 2 9 


■ ■ ■!( 







-0.06 -0.04 -0.02 



0.02 0.04 0.06 




^a jV (2 9 ), jv = 


--2 1 


+ u N (2<>),N-- 


--2 g 


+ u N (2 9 ),N = 


--2» 


-^-Sjv(2 10 ), JV 


= 2 7 


— 5«(2 10 ), JV 


= 2 8 


— S W (2 10 ), JV 


= 2 9 


^S w (2"), JV 


= 2 7 


^u N (2 n ), JV 


= 2 8 


-^5jv(2 n ), JV 


= 2 9 


MM.U 





0.72 0.74 0.76 0.78 



Figure 4: Numerical solutions un,\ for different values of N and A. The figures are 
zoom-ins of the solutions peaks. Results are provided for the test cases lb (top) 
and lib (bottom). 



in ( 19 ) is seemingly sharp. 





la 


Ila 


la 


Ila 


A 




rate 




rate 


ll e IU 2 (n) 


rate 


ll e IU 2 (f2) 


rate 


2 3 


l.lle-02 




3.37c-03 




1.12c-02 




3.47e-03 




2 4 


3.90c-03 


1.51 


1.19c-03 


1.50 


3.92c-03 


1.51 


1.22e-03 


1.50 


2 b 


1.37e-03 


1.51 


4.19e-04 


1.50 


1.38e-03 


1.51 


4.32e-04 


1.50 


2 K 


4.83e-04 


1.51 


1.48e-04 


1.51 


4.87e-04 


1.51 


1.52e-04 


1.51 


2 V 


1.69c-04 


1.52 


5.16e-05 


1.51 


1.70c-04 


1.52 


5.32e-05 


1.51 



Table 4: Dependence on the interaction radius A of the "errors" |||e||| and ||e||£2(Q) 
and the rate of convergence of approximate nonlocal solutions. Results are provided 
for the test cases la and Ila. 
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Figure 5: Numerical solutions un for different values of h and A show the conver- 
gence to the solution u of the fractional Laplacian as N — » oo and A — » oo. The 
figures are zoom-ins of the solutions peaks. Results are provided for the test cases 
lb (top) and lib (bottom). 





lb 


lib 


lb 


lib 


A 




rate 




rate 


ll e IU 2 (fi) 


rate 


ll e IU 2 (n) 


rate 


2 3 


1.41e-01 




6.11e-02 




1.42c-01 




6.27e-02 




2 4 


7.58e-02 


0.89 


3.40e-02 


0.85 


7.64e-02 


0.90 


3.49e-02 


0.84 


2 b 


4.16c-02 


0.87 


1.91e-02 


0.83 


4.19e-02 


0.87 


1.95c-02 


0.84 


2 6 


2.29c-02 


0.86 


1.07c-02 


0.83 


2.30c-02 


0.86 


1.09c-02 


0.84 


2 V 


1.24e-02 


0.87 


6.00e-03 


0.83 


1.25e-02 


0.88 


6.06e-03 


0.85 



Table 5: Dependence on the interaction radius A of the "errors" |||e||| and ||e||x,2(fi) 
and the rate of convergence of approximate nonlocal solutions. Results are provided 
for the test cases lb and lib. 
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6 Concluding remarks 

In this paper, we study a nonlocal diffusion operator which has as a special case 
the fractional Laplacian operator (— A) s ; exploiting a nonlocal vector calculus, we 
show that the solutions of nonlocal problems for the operator C converge to the 
solutions of fractional Laplacian problems as the nonlocal interactions become infi- 
nite. Through several numerical examples, we compare the nonlocal solutions with 
the solution of the fractional Laplacian equation on bounded domains, illustrate the 
theoretical results, and show that solving nonlocal problems is a viable alternative 
to solving fractional differential problems which feature infinite- volume constraints. 

A possible extension of this work consists in improving the rate of convergence, 
with respect to the grid size, of finite element approximations. The analytic solution 
of the fractional Laplacian equation (see e.g. (42 1) presents an abrupt change at 
the end points of the domain ft; in fact, it has a singular derivative in a; = —1,1. In- 
corporating the form of the singularity in the numerical scheme would yield greatly 
enhanced convergence; among the possible approaches we mention the singular ba- 
sis function method [201130] . In this approach, a set of supplementary functions that 
reproduce the functional form, assumed known, of the solution singularity is added 
to the ordinary finite element basis functions. Given the asymptotic expansion of 
the analytic solution as x —* — 1, 1, u(x) = 2]fc=i a kX Xk , where and Afe are the 
singular coefficients and exponents respectively, the supplementary basis functions 
have the general form 

ip k (x) = b{x)x x \ 

b{x) being an optional blending function. 

Another possible follow-up of the present work is to treat the nonlocal time- 
dependent diffusion equation 



Ut 



ft 



«(-,<) = inR"\n 
it(-, 0) = wq in ft. 



(45) 



This is a problem of interest in applications involving jump processes |10) . In fact, 
(451 is the master equation for a jump process, i.e., the deterministic equation that 



determines the time evolution of the probability density function restricted to a 
bounded domain ft. As done for the steady case, the solution u of (45) can be 
approximated by the solution u of 



Ut = Cm in ft 
2(-,i) = G infix 
0) = M in ft, 



(46) 



where C is defined as in ^3] A fractional time derivative might be considered as 
well. 
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Although in this work fractional Laplacian problems posed on all of R™ are 
not considered, such problems are indeed of interest so that the role of nonlocal 
diffusion operators for such problems would be interesting to explore. In this case, 
one would want to study the behavior of solutions of such problems as the domain 
ft, which for computational purposes has to be chosen to be finite, increases in size. 

We have only considered continuous Galerkin discretizations of the nonlocal 
diffusion problem, e.g., we have used continuous piecewise- linear finite element 
bases. For s < 1/2, discontinuous Galerkin methods are also conforming, i.e., 
satisfy (25). Thus it would be of interest to study such discretizations, especially in 
view of the fact that for s 1/2, both nonlocal diffusion and fractional Laplacian 
problems admit solutions containing jump discontinuities. 

Even if preliminary in nature, our numerical results presented here suffice to 
illustrate the theory and show the viability of using volume-constrained nonlocal 
diffusion problems as a means for defining approximations of fractional Laplacian 
problems posed on finite domains. Obviously, a natural follow-up of this work 
is to extend the numerical simulations to two-dimensional and three-dimensional 
settings for which computations would be even more challenging. 
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